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Abstract 

We present the results of non-equilibrium molecular dynamics simulations for 

the growth of a solid binary alloy from its liquid phase. The regime of high 

pulling velocities, V , for which there is a progressive transition from solute 

segregation to solute trapping, is considered. In the segregation regime, we 

recover the exponential form of the concentration profile within the liquid 

phase. Solute trapping is shown to settle in progressively as V is increased 

and our results are in good agreement with the theoretical predictions of Aziz 

[J. Appl. Phys. 53, 1158 (1981)]. In addition, the fluid advection velocity 

is shown to remain directly proportional to V , even at the highest velocities 

considered here (y ~ lOms^^). 

PACS numbers: 81.30.Fb, 68.45.-v, 02.70.Ns 
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I. INTRODUCTION 



Directional solidification (DS) of binary alloys is a reference experimental method to 
conduct carefully controlled tests of industrial casting. Besides their practical interest, DS 
experiments also bring insight into the fundamental study of basic instability morphologies 
of solid-liquid interfaces, such as cells or dendrites [I|]. By the combination of theoretical 
and numerical methods, much progress has been now achieved toward solving this difficult 
physical problem. However, most of the theoretical effort has concentrated so far on con- 
tinuous descriptions of the underlying phenomena [^j. In contrast, attempts to attack the 
problem at the atomistic level remain very few. To this extent, the present study may be 
considered as a first step to bridge the gap between both points of view. 

Because of the micrometer size of the growth structures, a direct and quantitative atom- 
istic simulation is still not at hand on this scale. Nevertheless, in the limit of large pulling 
velocities, a microscopic technique like Molecular Dynamics (MD) can still be used to follow 
atomistic phenomena occuring close to the solid-liquid interface. On one hand, MD was 
used to simulate laser-pulsed melting, for which the velocity is not controlled, but rather 
governed by heat diffusion. On the other hand, a first MD simulation of directed 

growth has been recently reported by Coura and al. 0. They considered the growth of a 
solid from a fluid phase with a density ten times smaller, a case which is probably more 
relevant to deposition from a vapor. 

In this paper, we present non-equilibrium molecular dynamic simulations of directional 
solidiflcation in two dimensions. We restrict ourselves to the case of rapid solidiflcation with 
a large temperature gradient, for which the interface remains plane on the atomistic scale, 
so that microstructures will not be considered here. After the description of the simulation 
details in section II, we present in section III the results obtained for the segregation pro- 
flies, the segregation coefficient and the advection velocity. In the last section, we finally 
discuss how this nonequilibrium simulation can be extended in the future to study different 
microscopic mechanisms involved in DS. 
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II. SIMULATION DETAILS 



Two atoms i and j separated by a distance r interact via the well-known Lennard- Jones 
potential : 

u^Ar) = Ae^AifY' - (1) 

The interaction energies between pairs of solvent-solvent, solute-solute and solvent-solute 
atoms are respectively denoted by Un, U22 and U12. The solvent potential parameters are 
chosen in order to describe Argon properties (en = 120K, an = 3.405A). For the solute, we 
take €22 = O.Seii and (722 = ctii- The cross-species parameters are fixed using the Lorentz- 
Berthelot rules : 

ei2 = \/^11^22 (2) 

and 

C^12 = + (^22)/2. (3) 

All the interactions are truncated at a cut-off radius = 2.5crii and the equations of 
motion are integrated using the Beeman algorithm [^] with a time step 6t = O.Olps. The 
particles coordinates are defined in a reference frame moving at the pulling velocity V in the 
X direction and periodic boundary conditions (PEG) are applied in the x and y directions. 
After integration of the dynamical equations over a time St, pulling is implemented by adding 
an increment —V6t to the x coordinate of each atom. In the reference frame, the solid-liquid 
interface is thus immobile when the stationnary state is reached. 

To simulate heat transport from the furnace to the system, four regions of fixed temper- 
ature are used (Fig. 0). Regions I, II, III and IV are centered at fixed positions, xj, xji, 
xjii and xiY and have a width of 20A for regions I and IV and lOA for regions II and III. 
In each region the temperature is kept constant by using a classical velocity rescaling. To 
maintain the solidification front between region II and III, we impose Ti = Tu < Tm and 
Tjjj = Tjv > Tm, Tm being the melting temperature of the alloy. 
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In Fig. |l| we represent the temperature gradient obtained after equilibration in the 
simulation box. The large difference between Tjj and Tjjj and hence the large gradient 
permits to localize the interface easily. With such a high gradient, instabilities cannot 
develop, which ensures the stability of planar interfaces. For the same reason we also take 
a small width Ly for the simulation box in order to reduce the natural roughness of the 
interface. More realistic systems, as compared to experiments, would correspond to xj = xn 
and xiii = xiv, together with a larger value of Ly. Because of the PBC in the x direction, we 
simultaneously have a solidification and a melting front, as in a melting zone experiment. We 
concentrate here on the solidification front and we use a wide liquid zone [xm < x < xjv) 
to allow solute diffusion. Fixing the temperature in four regions instead of two reduces 
considerably the rescaling of velocities within the liquid region. This is helpful to suppress 
pertubations and artefacts during the computation of microscopic quantities. The density 
difference between the solid and the liquid (~ 20%) is sufficient to induce advection of the 
liquid towards the front. A part of the associated momentum is then transmitted to the 
solid layer which in turn acquires a translationnal motion in the x direction. To avoid this 
finite-size effect, we rescale to zero the mean velocity within the deeper part of the solid 
(region I). 

Our system contains about 2000 atoms, 10 percent of which are solute atoms. Its size. 
La, = 400A in length and Ly = 60A in width, is relatively modest as compared to the size of 
systems currently used in MD. The reason is that the computational effort is here essentially 
spent in the time length of the simulation. To obtain good statistics, the simulation has to 
be long enough to allow each atom to perform several solidification-melting cycles. Since ten 
cycles require a time of lOL^/V, for the slowest velocity studied here (V^ = lOcms"^), this 
represents a simulation time of 4 x 10~^s (4 x 10^ MD steps). To increase the performance 
of our code we then adopt the 'cell lists' method ^ in which the box is divided into cells 
with a size slighty larger than the cutoff radius Tc- 
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III. RESULTS 



A. Concentration profiles 

The partition ratio, or segregation coefficient, is defined as 

k = c:/cj, (4) 

c\ and c* being the concentrations at the interface, respectively in the hquid and the solid. 



With our choice of the potential parameters (ei2 < en), is expected to be less than one [jTO 
so that solute accumulates in the liquid near the front. We ffist illustrate how the simulation 
reproduces solute rejection. The pulling velocity is fixed to = lms~^ and the solute atoms 
are initially placed at random in the simulation box. Fig. ^ is a snapshot of the system in 
its initial configuration. After a time At = 5ns, which corresponds to a spatial translation 
of the furnace of approximatly -L^j/S, the second snapshot (Fig ^) shows that almost all the 
solvent atoms initially close to the front are incorporated into the solid. Conversely, because 
of their poorer solubility in the solid phase, a majority of the solute atoms remain in the 
liquid phase. As compared to a solvent atom, it takes a longer time for a solute atom to 
cross the interface and hence the solute concentration is higher at the interface, as expected 
here. To quantify this segregation, we compute the proffie concentration by averaging over 
10^ uncorrelated spatial configurations (Fig. |^). The diffusion equation for the solute in the 
moving frame reads 

— D ^'^^ +V^^ (5) 
dt dx"^ dx 

In the stationary regime, it gives the well known theoretical profile for the solute concentra- 
tion in the liquid, 

c{x) = c\ + (cj - c\) exp(-x//,) (6) 

where Is = D jV ys, the diffusion length and the origin of the x axis is placed at the front. 
A good agreement is found between this expression and our simulations, the best fit giving 



a diffusion coefficient Dfu = O.sA^ps"^. To verify tliis value, we independently compute 
the diffusion coefficient profile D{x), with the help of Einstein's relation between the mean 
square displacement of the atoms and D, 

< \r{0) - r{t)\^ >= ADt. (7) 

In Fig. ^, we see that the diffusion coefficient increases abruptly between zero in the solid 
phase to roughly 0.4A^ps~^ in the liquid near the interface. This variation of D takes place 
over a distance of 30A that can be considered as a good aproximation of the interface width. 
In Eq. (5), D was assumed constant. We can check that the fitted value Dfn approximately 
corresponds to the average of D{x) within the liquid region near the interface. The same 
analysis is repeated for different pulling velocities between V = O.lms"^ and V = 5ms~^. 
For each velocity, we compute the concentration profile and extract Dfu from the diffusion 
length. In Fig. ||, Dfu is plotted as a function of the pulling velocity. As discussed below, 
the segregation coefficient tends to unity for the largest V and the calculation of Dfn is thus 
restricted to velocities V < 3ms^^. A reasonable agreement is found in each case between 
Dfit and the average value of the D{x) profile. 

In Fig ^, we plot the two concentration profiles obtained for V = lms~^ and V = 0.5ms~^. 
As expected, the smaller the velocity, the larger Is. If the agreement with the exponential law 
is good for the largest velocity, we notice that it is less satisfactory for the smaller one. This 
can be understood by a second look at Fig. ^ We observe that the thermal gradient induces 
an increase of the diffusion coefficient, as we go deeper in the liquid. As a consequence, we 
also have an increase of the diffusion length, which explains the disagreement observed for 
the largest x values. 

B. Segregation coefficient 

Now that we verified the ability of our method to reproduce both the segregation at the 
interface and the dependence of the diffusion length on the pulling velocity, we examine the 
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behavior of the segregation coefficient as a function of V. We see in Fig. ^ that the hquid 
concentration at the interface, cj, is larger for a lower value of V. This can be easily explained 
by the fact that, at large pulling velocities, the solute does not have enough time to fully 
diffuse over the interface width before it is incorporated into the solid. As a consequence, 
some solute is trapped in the solid and segregation at the interface is lowered. Several models 
have been proposed to describe solute-trapping HTll-p!^. The usual qualitative criterion for 
segregation is to compare two different characteristic times. One, ty = Si/V, corresponds to 
the growth of the solid over a distance comparable to the interface width. Si. The second, 
to = Sf/D, is the time needed for a solute atom to diffuse over the same distance 6i. The 
Aziz model predicts then 



'-T^ (8) 

for the segregation coefficient. Here kg is the segregation coefficient at equilibrium (i.e., 
when V = 0) and 

/3 = to/tv = 5i/ls (9) 

is the ratio of the two characteristic times (or lengths) defined above. Two remarks arise 
at this point. The ffist one concerns the neighborhood of the absolute stability threshold, 
V = Va, where the diffusion length becomes comparable to the capillary length, 1^ ~ d^. 
When (io ^ ^ii one also has /3 -C 1, so that ~ /^g. In this case, the absolute stability 
threshold is well described by the usual linear stability analysis (see [0). This is not true 
anymore, when do becomes comparable to as it is very likely the case for some reference 
alloys |T^. The second remark is that the assumption of continuous growth, on which 



Eq. (8) relies, is actually verified by the rough solid-liquid interfaces produced with our 
Lennard- Jones potential. 

In Fig 1^ we plot our estimates for A; as a function of the pulling velocity. The line 
is a best fit to the Aziz model. The parameter (3 is estimated by using 5i = 30A, the 
characteristic width of the interface extracted from the D{x) profile in Fig.^ and the value 



of the diffusion coefficient, D = 0.2A^ps~^ is estimated at the center of this interface region. 
The agreement between the simulated values and the model is rather good for a fit with a 
single free parameter. We obtain the value = 0.503 which is difficult to confirm because 
it is not possible to simulate segregation at much lower velocities with current computing 
power. Simulations at a fixed pulling velocity thus provide an indirect way of determining 
fcg. The fact that the interface is moving permits all atoms, and especially the solute ones, 
to be part of the fluid phase and then to speed up relaxation to the equilibrium situation. 
For an immobile interface, one would have to wait for the very slow diffusion of solute atoms 
within the solid phase. The solid-liquid equilibrium for 3D binary Lennard- Jones mixtures 



has been determined by Monte-Carlo simulations [|T0],|T5|. It would be interesting to have 
similar 2D results available to test our indirect evaluation of k^. 



C. Advection velocity 

We finally study the advection of the fluid near the solidification front. The Lennard- 
Jones mixture considered in this work produces a substantial density difference between the 
solid and the liquid. This is illustrated in Fig. |], where we plot the density profile along the 
X direction perpendicular to the interface. From these data, we estimate both densities at 
the interface, pi and p]. We find that they do not depend on the pulling velocity and that 
the normalized density difference is also a constant, {pi — p])/ p] — 0.24. This rather high 
value for a Lennard- Jones potential is due to the large temperature gradient imposed at the 
interface. As a consequence, the solid cannot grow unless the liquid is advected toward the 
front with a velocity Vad- In the laboratory frame, the solid is at rest, while the liquid has 
an advection velocity Vad all along the x axix (Fig. Mass conservation at the interface 
imposes that V^^ varies with V as 

= ^^H^y- (10) 
Pi 

A linear variation is effectively recovered in our simulations (Fig. |r^), with a slope comparing 
well to the estimate extracted from the density profile (see Fig. H). It is interesting to note 
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that there is no deviation from the hnear law, even at the highest velocities studied in 
this work. On one side, one could have expected that the onset of solute trapping would 
modify the densities. This is not the case here, probably because the solvent and solute 
atoms have the same atomic radius, (J22 = <^ii- On the other side, a deviation would also 
appear as V is increased, if the solid was created with more and more defects. We would 
then have a decrease of V^^ because the density of the solid would also decrease. This type 
of deviation is not observed any more in our simulations. We finally see in Fig. ^ that 
the advection velocity increases in magnitude with x. As we go deeper in the liquid, the 
density progressively decreases in the thermal gradient and an extra advection velocity is 
thus necessary to compensate the density drop between two adjacent slices of alloy. 



IV. CONCLUSION 

With the help of MD simulations, we directly observed and analyzed quantitatively two 
important phenomena occuring at the atomistic level in rapid DS. The first one is the cross- 
over between the regime of solute segregation and the regime of solute trapping, which 
sets in as the pulling velocity, V, is increased. The good quantitative agreement of our 
estimates for the segregation coefficient, k{V), with existing theories [irTHT^ confirmed that 



MD is a valuable simulation tool for DS. The second point is the advection of the liquid 
phase towards the solidification front. The advection velocity at the interface was shown 
to remain proportional to V, in the whole range of values considered here. This effect is 
more likely masked by stronger phenoma like convection in the usual experimental set-ups. 
It may however become relevant in micro-gravity or for thin sample DS. 

The use of Lennard- Jones potentials is rather restrictive to compare the present results 
with experiments. In the future, potentials based on empirical descriptions, like the embed- 
ded atom method |T^, the glue model or the effective medium theory |jl8[| will be used 



to study technical materials like metals. This will allow us to explore important issues like 
anisotropy and facetting, or the role of defects and constraints. In parallel, it will be neces- 
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sary to simulate more extended systems, in order to explore the front instabilities from an 
atomistic point of view. In this last category, we can cite the oscillations of the solidification 
front just above the absolute stability velocity pJ|j20[] . 
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FIGURES 

FIG. 1. Temperature profile along the x axis. The different regions of the simulation box in 

which the temperature is kept fixed are bounded by the dotted lines. The melting temperature is 

FIG. 2. Snapshots of the system at two different times. Large and small radius respectively 
correspond to solvent and solute atoms, a) At t = 0, solute atoms are randomly placed in the 
simulation box. One slice of liquid in contact with the interface is marked (dark grey), b) At time 
t = 5ns, almost all the solvent atoms have solidified while a majority of the solute atoms remained 
liquid. 



FIG. 3. Concentration profile c{x) along the pulling direction, perpendicular to the interface. 

FIG. 4. Diffusion profile D(x) along the direction perpendicular to the interface. The horizontal 
dotted line is the average diffusion coefficient obtained by fitting the concentration profile of Fig. 
3 to an exponential law. 



FIG. 5. The average diffusion coefficient as a function of the pulling velocity. 

FIG. 6. Solute concentration profiles obtained for V = lms~^ (circles) and V = 0.5ms~^ 
(diamonds). The continuous lines are best fits to Eq. (6). 

FIG. 7. Variations of the segration coefficient with the pulling velocity. The solid line is the 
best fit to the Aziz equation. 

FIG. 8. Density profile p{x) along the direction perpendicular to the interface. The upper and 
lower arrows indicate respectively pi and p| at the solid-liquid interface. 

FIG. 9. Velocity profiles, Vx{x), along the direction perpendicular to the interface. Open and 
filled circles respectively correspond to V = 3ms~^ and V = Sms"^ 
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FIG. 10. Interface advection velocity V^^, plotted as a function of the pulling velocity V. The 
solid line corresponds to Eq. (10). 
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